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The addition of suitable volume forces to the Navier-Stokes equation allows to simulate flows in 
the presence of a homogeneous shear. Because of the explicit form of the driving the flows are 
accessible to rigorous mathematical treatment and to accurate quantitative modelling in their 
global properties. The statistics of the fluctuations provide insight into generic behaviour of 
non-equilibrium systems and into the presistence of anisotropies at small scales. Correlation 
functions can be used to identify dominant large scale dynamical processes that are relevant 
for most of the momentum transport across the shear. The numerical studies of homogeneous 
shear flows complement analytical and experimental investigations and contribute to bridging 
the gap between ideal homogeneous isotropic turbulence and the more realistic heterogeneous 
turbulence. 

1 Introduction 

The turbulent flow of a fluid is a paradigmatic example for a non-equilibrium, nonlinear, 
and self-organizing system. The flow has to be driven continuously to compensate the 
viscous friction due to shearing; the advection of the flow field by itself introduces a char- 
acteristic nonlinearity; and the flow shows a hierarchical organization of smaller whirls on 
top of larger ones 1 . It is one of the aims of turbulence theory to derive from the equations 
of motion the laws that govern this self-organization in space and time. Partial progress has 
been made for the theoretically most appealing case of homogeneous, isotropic turbulence, 
where we have the Karman-Horwath equation for third order moments 1 , a mean field the- 
ory 2 , and a variety of models 13 , but for the most part an analytical theory is not within 
reach. In the absence of such a formalism numerical simulations provide both access to 
the properties of the equations as well as guidance in developing intuition. 

The power laws that describe the distribution of energy on the various scales quite 
well can be obtained by dimensional arguments, following Kolmogorov, Onsager, v. 
Weizsacker, Heisenberg and Richardson 1 . These power laws are truncated on large scales 
by the forces that stir the fluid and on small scales by the onset of viscous dissipation. 
Therefore, one would like to study huge systems with infinitesimal viscosity in order to 
obtain as large a scaling range as possible. In terms of the large scale Reynolds number 

Re=UL/v (1) 

formed by an external velocity scale U, an external length scale L and the viscosity of the 
fluid v, this is the limit of Re — > oo. Phenomenological arguments show that the ratio 
between the externally imposed length L and the large eddy turnover time T = L/U to the 
smallest scales rj and r that appear in the dynamics scale with Reynolds number like 3 

r)/L - Re- 3 / 4 t/T ~ Re" 1 ' 2 . (2) 

Thus, the number of modes that are needed in order to represent a 3-d flow field increases 
like Re 9 / 4 . Since the maximal time step in the integration is not determined by the time 
scale of the flow, but by stability considerations following from the spatial discretization, 



the number of time steps increases like Re 3 / 4 , so that the total operation count increases 
like Re 3 \ As a consequence, the largest simulations run to date, on the Earth simulator 
with 4096 3 modes, achieve a Reynolds number of about 230.000 only, and follow the flow 
field for just one large scale time unit 4 . However, the comparison of data by Sreenivasan 5 
suggests that already with a resolution of about 256 3 or 512 3 one can reach into the iner- 
tial range behaviour and thus begin to study the dominant statistical behaviour of turbulent 
flows. This then opens up the possibility to study flows more realistic but also more com- 
plicated than homogeneous isotropic turbulence. 

In contrast to the theorists ideal, realistic flows are driven by anisotropic forcings (e.g. 
a uni-directional winds), influenced by boundary layers (near walls) and are in many other 
respects far from the case of homogeneous, isotropic turbulence. A flow that is interme- 
diate between realisitic flows and the ideal cases, and that can be used to bridge the gap 
is turbulence with a homogeneous, linear shear profile superimposed. It mimicks the situ- 
ation that one expects to find, at least locally. An investigation of such flows then allows 
to study the interaction between the turbulence and the shear, the range over which the 
anisotropics affect the scaling behaviour, and the consequences of the shear for dynamical 
transport processes in the flow. 

The problem with homogeneous shear flows is that it is not obvious how to maintain 
the shear. Driving the fluid from the boundaries, e.g. by oppositely moving side walls, will 
produce the familiar boundary layers with large gradients and a bulk region with reduced 
ones. Rogallo 6 thus suggested to drive the fluid by continuously shearing the computational 
grid. Obviously, after a while, the grid will be extremely distorted and a remeshing will 
be necessary in order to restore the initial resolution. When the appropriate terms are 
introduced into the Navier-Stokes equation, this amounts to a periodic driving of the fluid, 
and thus a temporally inhomogeneous situation. 

In 1999 we started a research programme based on a novel method with which we can 
simulate homogeneous shear flows with prescribed forcing, shear, and residual turbulence. 
The essential idea is to employ a body force that can be tuned to fit the appropriate re- 
quirements. It is a numerical tool, not realizable in experiments, that is fairly versatile and 
amenable to mathematical analysis and provides insights impossible to obtain otherwise. 

The outline of the article is as follows. We begin in section 2 with a description of 
the numerical methods used. We then turn in section 3 to a description of the global 
properties of the flow, the connection to rigorous mathematical bounds and to a model 
for the relaxation of the turbulent energy. The fluctuations around the mean properties, 
their probability distributions and their scaling with Reynolds number are discussed in 
section 4. The dynamical properties as reflected in dynamical correlation functions and 
their significance for large scale momentum transport are discussed in section 5. A brief 
outlook concludes the article. 



2 Simulating homogeneous shear flows 

The Navier-Stokes equation for an incompressible fluid driven by a divergence free volume 
force F is, in dimensionless form, 

-^ + (u-V)u= -Vp+— V 2 u + F, (3) 
V • u = , (4) 
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with p(x, t) the pressure, and u(x, t) the velocity field that satisfies V • u = 0. The difficult 
and time-consuming part of the integration is the evaluation of the nonlinear advection term 
(u • V)u. High Reynolds number simulations therefore prefer Fourier representations for 
the velocity field, so that Fast Fourier Algorithms can be used to obtain the gradients as a 
local multiplication in wave number space. The aliasing problem is handled with the 2/3 
rule 7 . In order to break the inherent periodicity of the Fourier modes, which implies that 
to every region with positive shear there will be one with negative shear shifted by half 
a period, we bound the domain by free-slip walls in the shear direction: these boundary 
conditions are then compatible with a Fourier representation of the velocity field. 

The coordinates are chosen such that x points downstream, y in the direction of the 
shear, and z in spanwise direction. In x and z we take periodic boundary conditions, and 
in the y direction free-slip conditions, i.e., 

, du x du z 

u y = 0, and = — = (5) 

at top and bottom surfaces. The size of the domain is [0, 2-k] x [0, L y ] x [0, 2ir]. The length 
L y is taken to be 1 or tt. The typical Fourier resolution is 256 x 129 x 256. As a rule of 
thumb the spectral resolution is considered sufficient when the largest wave number k max 
and the smallest turbulent scales i] satisfy k max i] > 1, although for quantities sensitive 
to gradients the fluctuations in r\ may have to be taken into account, so that even more 
restrictive conditions have to be satisfied (Schumacher and Sreenivasan, work in progress). 
Because of the dealiasing the largest wave number is given by V2N/3 for N modes 7 . 

For turbulent flows it has become customary to use not the external Reynolds number 
Re but an easily computable and measurable intrinsic Reynolds number R\ = U rms \jv. 
It is based on the root mean square velocity U rms and the Taylor length scale A obtained 
from the gradients of the velocity field, 

2 ((du x /dxf) 

= <«!> ' (6) 

In isotropic turbulence the relation between the large scale Reynolds number Re and the 
Taylor-Reynolds number R\ is given by 3 

ReKQ.lbR\. (7) 

Our simulations reach up to R\ w 150 on a 512 3 grid for the bounded flow domains, the 
earth simulator 4 manages R\ w 1200 for isotropic turbulence. 

The flow is naturally homogeneous in x and z, so that for statistical purposes we can 
form averages over planes y — const, henceforth denoted (. . .)a- It cannot be fully homo- 
geneous in the y-direction because of the vanishing gradients at top and bottom surfaces, 
but the simulations show that the widths of these free slip boundary layers is small and 
decreases with increasing Reynolds number (see the inset of Fig. 1 and Ref. 7). 

Two forms of driving are typically used: we can fix a force field F and study the flow 
that results, or we can prescribe the mean profile and adjust the force so that this mean 
profile is maintained. In the latter case, the force compensates the contributions that would 
result from the nonlinear term and the pressure, in order to keep the amplitudes of a set of 
preselected Fourier components constant. 

The mean profiles and the statistical properties of this flow were compared to exper- 
iment and other simulations in 9 . It turns out that the method easily yields a statistically 
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stationary state and that the dynamics of the flow is much less violent than in the case of 
the Rogallo remeshing, i.e. the bursts in energy and enstrophy typically found there 10 ' 11 
are absent. 



3 Globally averaged properties 



A quantity that is amenable to rigorous mathematics without supplementary hypothesis 
and additional models is the total dissipation in the mean, i.e. the turbulent dissipation 
averaged over the flow domain and over time, 



u(\Wu\ 



>V,T ■ 



(8) 



The optimal bound theories of Busse 12 , the variational approaches of Constantin and Do- 
ering 13 , Kerswell 14 and Nicodemus et al 15 and their recent extension to volume driven 
flows 16 ' 17 show that this dissipation is rigorously bounded from above by an expression of 
the form 



e < civ- 



il 2 



+ c 2 - 



(9) 



e • ~* t 

with numerical coefficients c\ and c 2 . The velocity scale U rms is set by the root mean 
square of the velocity and I is the length scale of the external driving. If a Reynolds 
number is formed with U rms , i.e. Re = U rms £/v, then this bound is consistent with the 
expectation that e ~ Re 3 for large Re. 

For the specific case of a sinusoidal driving the bounds and the results from a numerical 
simulation are shown in Fig. 1. The coefficient from the upper bound is c 2 = 7r 2 /V216 « 
0.67, whereas the numerical simulations indicate c 2 ~ 0.2, about a factor 3 lower. It is 
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Figure 1. Dimensionless dissipation ratio, /3, vs. Reynolds number for a flow driven by a constant shear force. 
The results of the direct numerical simulations are indicated by diamonds with error bars from the statistics of the 
fluctuations. The straight dotted line indicates lower limit to the dissipation, obtained for a laminar shear flow. 
The three other lines are analytically derived and succesively improved bounds. The inset shows the mean flow 
profiles (u x )a,t vs y for two different Reynolds numbers. 
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interesting to note that the difference between the bounds and the numerical simulations is 
much smaller than in the case of wall driven shear flows 15 , where the difference is about 
a factor 10, and that in contrast to that case the dissipation does not tend to decrease with 
increasing Reynolds number. This suggests that volume forces provide a more efficient 
stirring than moving boundaries. 

The instantaneous dissipation e(i), energy content E(t) = (u 2 )y/2, and energy uptake 
I(t) = (u • F)v are not constant but fluctuate fairly irregularly. An equation for their 
dynamics can be derived with less mathematical rigor but still directly from the Navier- 
Stokes equations with minimal assumptions. Multiply the Navier-Stokes equation (3) once 
with u and once with F and average over the fluid volume. Then: 



f = -e(t) + /(t), 
g = -(F.[(u-V)u])v 



Wu-AFW + (F 2 



>v ■ 



(10) 



(11) 



In case the force field is an eigenfunction of the Laplacian, the term (u ■ AF) y becomes 
proportional to I(t). The equations can then be closed by relating the energy dissipation 
and the term quadratic in the velocity in the second equation to the energy content, e = 
c d E 3 / 2 , and (F ■ [(u ■ V)u])y = CfE. Then 

AE 

dl 
dt 



-c d E 3 ' 2 (t)+I{t) 
-c f E(t) - vXI(t) 



F. 



(12) 
(13) 



where the last term contains the norm of the force profile, F = (F 2 )y. The model con- 
stants Cd and Cf can be determined from turbulent flow simulations. These equations have 
a stationary state, corresponding to the time average energy content and energy uptake. 
Interestingly, the relaxation to this stationary state is oscillatory. Fig. 2 shows a typical 
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Figure 2. Comparison of the energy model with direct numerical simulations of a shear flow. The left panels 
show time traces of the energy content (top) and energy uptake (bottom). The panel on the right shows a typical 
trajectory from the DNS, superimposed on the color contours of the probability density from the model in the 
presence of additive white noise. 



5 



trajectory from a numerical simulation: the fluctuations that are superimposed on the mean 
values always drive the system away from the stationary state, but the relaxation towards 
the stationary state is not unlike the dynamics in the model. When the fluctuations are 
modelled as additive white noise the colored probability density results. The intensity of 
the noise was determined from the requirement that the second moments of the measured 
fluctuations be reproduced. 

The model provides a promising starting point for the analysis of relaxational behaviour 
in excited turbulent flows (previously observed in 18,19 ) and for the response dynamics of 
periodically driven flows 20 ' 21 . In view of the exact results that could obtained for the 
stationary situation it will also be interesting to see whether similar results can be obtained 
for periodically driven flows. 

4 Local fluctuations 

Simulations and model in the previous section show that energy dissipation and energy 
uptake in a turbulent flow are fluctuating quantities, and even though they are equal in the 
mean, they can be different for short times. The distribution of such fluctuations is a matter 
of debate, with experiments by Pinton et al. 22 showing strongly non-Gaussian fluctuations, 
whereas the data of Cadot et al. 23 are compatible with a Gaussian distribution. Since 
difference may be related to the size of the domains over which the flow field was averaged 
we studied the local energy uptake. This also shows clear signatures for the presence of 
negative energy uptake, indicating a reverse flow in energy from the turbulent fluid to the 
stirrer. Such events are very rare, and if the average over the fluid volume is taken, they 
essentially never happen. However, locally the fluctuations of energy uptake rates can 




Figure 3. Statistics of energy uptake for different vertical positions in the shear layer. The probability density 
functions (PDF) normalized by the mean values are shown in the left panel. They are based on more than 100 
turbulent snapshots separated in time by 1.5 large eddy turnover times with a total of about 6 X 10 8 data points. 
The PDF's are shown for three different Reynolds numbers, Re = 270 (red), 360 (blue) and 390 (green) and 
two different positions yo in the normal direction (thin lines for 0.27 and thick lines for 0.02). The right panel 
shows the ratio of the probability for positive and negative energy dissipations. 
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take on both signs, and quite frequently become negative. Statistics of the instantaneous, 
pointwise energy injection rate e in (x, t) = u(x, i)-F(x) were analyzed in 24 . Since the 
system is invariant under translation in downstream and spanwise direction, but not in the 
normal direction, we study the distributions for planes parallel to the bounding free-slip 
surfaces separately. The probability density functions of the energy input rate in units of 
its ensemble average, (e in ) v . T , and for different positions between the plates are shown in 
Fig. 3. 

The probability density functions for different Reynolds numbers collapse nicely when 
normalized by the mean, but the distributions vary considerably across the layer. Negative 
values, i.e. transfer of energy to the stirrer, occur fairly frequently, and are more likely 
further away from the center. The relative frequency between positive and negative energy 
uptake shows an exponential relation. Such exponential relations have recently been found 
in many non-equilibrium systems where they could be connected to fundamental symmetry 
properties of non-equilibrium invariant measures 25 . For hydrodynamic systems these ideas 
do not strictly apply, since the Navier-Stokes equation is not reversible. But the study of 
Farago 26 shows that even in the absence of that symmetry, as e.g. for a Brownian particle, 
similar relations can be justified analytically. 

Since homogeneous shear flows are in a sense the first step away from ideal isotropic 
turbulence, one can ask how the statistical properties at the smallest scales of the turbulent 
flow are affected by the shear. In a pioneering paper Lumley 27 predicted a rapid R^ 1 decay 
of such anisotropics with Reynolds number. Recent systematic measurements in simple 
shear flows 28-31 for Taylor microscale Reynolds numbers up to R\ <~ 10 3 show clear 
deviations from the predicted decay. Moreover, direct numerical simulations (DNS) for 
moderate Reynolds numbers confirm a persistent anisotropy and reveal a relation to typical 
large-scale flow structures, so-called streamwise streaks 10 ' 9 ' 32 ' 33 . While the anisotropics 
still tend to become smaller with increasing Reynolds number, the decay is slower than 
R~ x \ 

5 Dynamical correlations 

Various kinds of structures in shear flows such as vortices, streaks or waves have been 
identified and considerable efforts have gone into identifying their dynamical origins and 
evolution. The isosurfaces of the downstream component of the flow field in Fig. 4 show 
that even though the driving is homogeneous in spanwise direction the flow organizes 
into large regions with very high positive (red) or negative (blue) downstream velocity. 
Such structures are called streaks. They induce strong gradients and shear instabilities 
and play an important role in a turbulent recycling process described by Waleffe 34 that 
consists of three steps: i) downstream vortices mix fluid in the normal direction and drive 
modulations in the downstream velocity, forming streaks; ii) streaks undergo an instability 
to the formation of vortices pointing in the normal direction; iii) the mean shear profile 
now turns these vortices again in downstream direction, thus closing the loop. Of these 
processes the ones in step iii) and ii) are reasonably fast, whereas the one in i) is fairly 
slow, since it is connected with a non-normal amplification and thus only linear in time. 

The indicator for non-normal amplification that we focus on here is a temporal cross- 
correlation function 35 . Since the vortex drives the streak a cross-correlation between the 
vortex and the streak should be asymmetric in time: if the streak is probed after the vortex 
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Figure 4. Isosurface plot of a turbulent streamwise velocity snapshot for a DNS of homogeneous shear turbulence 
at R\ = 170. Level at v x = 2 is coloured red, level at v x = — 2 blue. The elongated structures can be identified 
as streaks. 

then there might be a correlation, if it is probed before then there should not be a correla- 
tion. For the fluctuations in a linearization around a linear shear profile this can be analyzed 
analytically 35 . 

In order to see this in turbulent flows we study Eulerian spatial and temporal cross- 
correlation functions between the downstream and normal velocity components at fixed 
heights yo, i.e., 

C xy (Ax, At; y ) = (v y (x, y ,z, t)v x (x + Ax, y , z, t + At)) x<Zit , (14) 

The correlation function has elongated oval like isocontours that are not aligned with the 
axes (Fig. 5). The angle by which it is tilted gives a velocity that seems to differ a bit 
from the mean velocity of the flow (indicated by the green line). One-dimensional cuts, 
such as the one along the red line at fixed position, support an asymmetry, indicative of 
the non-normal amplification process 36 . The cut along the line given by the mean velocity 
(green) is what would be expected in a hot-wire anemometor at a fixed position using 
Taylors-hypothesis 37 . 

6 Outlook 

The investigations on homogeneous shear flows presented here have helped to understand 
the dynamics of energy and dissipation, of the behaviour of the fluctuations and the dy- 
namical process that are active. There are several directions in which one can proceed: 
one option is to add tracers in order to study the dynamics in advected Lagrangian frames, 
expecially the relative dynamics of three or more particles since it should contain informa- 
tion about the statistically conserved quantities in the cascade 38 , and since it connects to 
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Figure 5. Spatio-temporal cross-correlation function of a turbulent shear flow evaluated at a fixed normal position 
yo between the free-slip boundaries. The left panel shows the contour plots of C xy (Ax, At; j/o) as defined in 
(14). The right panel shows one-dimensional cuts in time for a fixed position (red) and along a position moving 
with the mean velocity, Ax = (u x )(yo)At (green). The green correlation function is one to be expected 
on the basis of Taylor's frozen flow hypothesis 37 . The inset shows the asymmetry parameter Q = (C(t) — 
C(-t))/(C(t) + C(-t)). 

the experiments of Bodenschatz et al 39 . Along similar lines, it should be interesting to add 
scalar fields in order to study their mixing dynamics 40 . Active particles, like long flexible 
polymers, will be stretched by the flow field 41 and will interact with the flow and will inter- 
fere with the large scale dynamics, so as to reduce turbulent drag in the fluid. The insights 
gained in the homogeneous shear flow studies will be extremely valuable in understanding 
the properties of these added fields. 
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